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Three-body continuum states are investigated with the hyperspherical method on a Lagrange 
mesh. The i?-matrix theory is used to treat the asymptotic behaviour of scattering wave functions. 
The formalism is developed for neutral as well as for charged systems. We point out some specificities 
of continuum states in the hyperspherical method. The collision matrix can be determined with a 
good accuracy by using propagation techniques. The method is applied to the ®He (=a+n+n) and 
6 Be (=a+p+p) systems, as well as to 14 Be (= 12 Be+n+n). For 6 He, we essentially recover results 
of the literature. Application to 14 Be suggests the existence of an excited 2 + state below threshold. 
The calculated B(E2) value should make this state observable with Coulomb excitation experiments. 

PACS numbers: 


I. INTRODUCTION 

Three-body systems present a large variety of interesting features Hi The discovery of a halo structure in 6 He 
Q triggered many experimental and theoretical works on exotic nuclei, such as 6 He, n Li or 14 Be. The bound-state 
spectroscopy of Borromean systems is now relatively well known. On the experimental side, current intensities of 
radioactive beams are high enough for precise measurements of spectroscopic properties, such as energies, r.m.s. radii 
on quadrupole moments. On the theoretical side, several methods have been developed, and provide accurate solutions 
of the three-body Schrodinger equation. 

The hyperspherical harmonic method (HHM) is known to be well adapted to three-body systems Hi. The six 
Jacobi coordinates are replaced by five angles, and a single dimensional coordinate, the hyperradius. The HHM trans¬ 
forms the three-body Schrodinger equation into a set of coupled differential equations depending on the hyperradius. 
It has been applied to many exotic nuclei. 

Recently, we have combined the HHM with the Lagrange-mesh technique [(J. The Lagrange-mesh method (see 
Ref. |j] and references therein) is an approximate variational calculation that resembles a mesh calculation. The 
matrix elements are calculated at the Gauss approximation associated with the mesh. They become very simple. In 
particular, the potential matrix elements are replaced by their values at the mesh points. In spite of its simplicity, 
the Lagrange-mesh method is as accurate as the corresponding variational calculation. This was shown for two-body 
Q as well as for three-body 0 systems. 

In the present work, we extend the formalism of Ref. 6] to three-body continuum states. The information provided 
by continuum states is a natural complement to the bound-state spectroscopy. Experimentally, three-body continuum 
states are investigated through breakup experiments (see for example Ref. 0 ). On the theoretical point of view, 
various methods have been developed. Some of them, such as the Complex Scaling Method Q, or the Analytic 
Continuation in the Coupling Constant jTf)l | deal with resonances only; they cannot be applied to non-resonant states. 
Other methods, such as the R-matrix theory 0 are more difficult to apply, but can be used for non-resonant, as well 
as for resonant states. 

Applications of the R-matrix method to two-body systems have been performed for many years in nuclear as well 
as in atomic physics. In nuclear physics, applications to three-body systems are more recent IT3 . The R-matrix 
theory allows the use of a variational basis to describe unbound states. It is based on an internal region, where the 
wave function is expanded over the basis, and on an external region, where the asymptotic behaviour can be used. In 
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three-body systems, the hyperspherical formalism is very efficient for bound states. For unbound states, however, it 
raises problems owing to the long range of the coupling potentials [l^ ■ In the f?-matrix framework this can be solved 
by using propagation techniques |~Hil |. 

In two-body systems, the Lagrange-mesh technique associated with the i?-matrix formalism has been applied in 
single- 0 and multi-channel [jjj calculations. The purpose of the present work is to extend the method to three-body 
systems. Another development concerns the application to charged systems. Many exotic nuclei are unbound, even 
in their ground states, due to the Coulomb force. We show applications to the a+n+n and 12 Be+n-|-n systems, for 
which two-body potentials are available in the literature. The mirror systems are also investigated. 

In Section 2, we summarize the three-body formalism, and present the i?-matrix method. Section 3 is devoted to 
applications to 6 He and 14 Be, with the mirror systems. Concluding remarks are given in Section 4. 


II. THREE-BODY CONTINUUM STATES 
A. Hamiltonian and wave functions 

Let us consider three particles with mass numbers A t (in units of the nucleon mass mjv), and space coordinates rq. 
A three-body Hamiltonian is given by 

3 3 

H = J2 Ti+ E V ^~ r ^ (!) 

i=l i>j =1 

where Tj is the kinetic energy of nucleon i, and Vij a nucleus-nucleus potential. We neglect three-body forces in this 
presentation. 

The HHM is known to be an efficient tool to deal with three-body systems. This formalism is well known, and we 
refer to Refs. il for detail. Startingfrom coordinates rq, one defines the Jacobi coordinates x k and y k (k = 1, 2,3). 
We adopt here the notations of Ref. [gj. The hyperradius p and hyperangle a k are then defined as 

P 2 = x k + yh 

a k =arctan—. ( 2 ) 

Xk 

The hyperangle a k and the orientations f l Xk and fl yfc provide a set of angles £l§ k . In this notation the kinetic energy 
reads 


3 

T p = Ti- T cm 
2=1 


( &_ 5 d K 2 (n 5k ) 

2rriN \dp 2 ^ p dp p 2 


( 3 ) 


In Eq. m,Tcm is the c.m. kinetic energy, and K 2 is a five-dimensional angular momentum [Tg1| whose eigenfunctions 
(with eigenvalues I\(K + 4)) are given by 

yK e L y M L (Os) = tfK* («) [Yt. (^) ® Y, y (H y )] LML , 

< /bE( a ) = {c°saY x (sinaY y pi, v+2 ’ ex+2 \cos2a), (4) 

where Pn a, ^\x) is a Jacobi polynomial and a normalization factor Q (here k is implied). In these definitions, 

K is the hypermomentum, (£ x ,£ y ) are the orbital momenta associated with ( x,y ), and n is a positive integer defined 

by 


n = {K - 4 - £ v )/ 2. 

Introducing the spin component x SMs yields the hyperspherical function with total spin J 

JM r a . 


yKL^)® x k 


( 5 ) 

( 6 ) 


where index 7 stands for (£ x ,£ y , L, S). 
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A wave function 
basis functions © as 


, solution of the Schrodinger equation associated with Hamiltonian ©, is expanded over 


7 K 


( 7 ) 


where n = (—1)^ is the parity of the three-body relative motion, and Xjk(p) are hyperradial wave functions which 
should be determined. Rigorously, the summation over (jK) should contain an infinite number of terms. In practice, 
this expansion is limited by a maximum K value, denoted as K max . For weakly bound states, it is well known that 
the convergence is rather slow, and that large K rnax values must be used. Typically 100 — 200 terms are necessary for 
realistic K max values. 

The radial functions X^k(p) are derived from a set of coupled differential equations 

' ^ (#_ _ C K (C K + 1 ) 

2m n \dp 2 p 2 

+ Vkj,k'j'(p) X^yk’(p) = 0 ) ( 8 ) 

K'y 



with Ck = K + 3/2. The potential terms are given by the contribution of the three nucleus-nucleus interactions 


Vk7,k'y(p) ~ Z(V£Z*(P) + Vk^k*y(p))> 


(9) 


i= 1 


where we have explicitly written the nuclear ( N ) and Coulomb (C) terms. 

Assuming the use of (* 1 , y^) for the coordinate set, the contribution i = 1 is directly determined from 


v£$M=fy$rmv 2 a (^) 


( 10 ) 


where = AiAj/(Ai+Aj). The terms * = 2,3 are computed in the same way, with an additional transformation using 
the Raynal-Revai coefficients ^j. Definition m is common to the nuclear and Coulomb contributions. Integrations 
over Da, and fl y are performed analytically, whereas integration over the hyperangle a is treated numerically. For the 
Coulomb potential, the p dependence is trivial; we have 


V vi^L(o) = zii. 


(ii) 


where , K , is an effective charge, independent of p, and calculated numerically from Eq. 1101 ) and from Raynal- 
Revai coefficients E3- Examples of matrices z J ' K are given in Ref. m for the cc+p+p system. Knowing the analytical 
p-dependence of the potential is crucial for continuum states (see below). Notice that, to derive Eq. (till) , one assumes 
the 1 /1r*j — r^| dependence of the Coulomb potential. Using a point-sphere definition is straightforward, as the 
difference can be included in the nuclear part. 


B. Asymptotic behaviour of the potential 

For small p values the potential must be determined by numerical integration of Eq. E3>. However, analytical 
approximations can be derived for large p values. For the Coulomb interaction, definition GH> is always valid. Let us 
now consider the nuclear contribution. After integration over D x and D y , a matrix element between basis states © 
is written as 
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T J-vt-vt’Jy ( \ r r C jtxtyf U/ {P COS a\ J'J 2 2 , 

V KL,K'L'\P ) = OLL'(>t y i' y J (p K ^a)V N y-j= j (p£,{a) sin acos acta 


. i rP (n . 1 

*i VK rCl~~ ~ 1 ' 


= K x y K* y 5 LL ,S iti e'^ 


(ty + o O 1 


xP ? T2 


“ s P 3 7o 


2—7 — 1 

p 


Cw + 2 >U + 2 ) 


1 - — 
P 2 




2 dw 


( 12 ) 


To deal with the spin, the coupling order in Eq. © is modified in order to introduce the total spin of the interacting 
particles j x — t x + S. This is achieved with standard angular-momentum algebra, involving 6 j coefficients. If the 
tensor force is not included, we also have £ x = i' x . For large p values, and if the potential goes to zero faster than 
1/m 2 , we can use the following expansions |jj) 


Pi“^( 2 x- 1 )=^ 




m =0 




(_ 1 )n+m T ((3 + n+l)T(a + (3 + n + m + l) 
m\{n — m)\ T(0 + m + l)T(a + /3 + n + 1) 


( 1 -*)“=£ 


m —0 


(-xy 


and we end up with the asymptotic expansion of the potential 


(13) 


yW^(p)« Sll'Sui. 


y t v 0 e x +e' x +3 0 . 

y k—0 H 


E 


Vk 

2 k 7 


where 


v k = M^ V K& y 


u tx+P x +2k+2y 


(- 2 =) 

\\/P23J 


du 


EM) 


k—m i —m 2 


P _L I 
l 2 / r 2 

k — 777-1 — m 2 


(A/ + 2 >^*+ 2 ) (^y+ 2 ’^x+ 2 ) 


(14) 


(15) 


Owing to the finite range of the potential, the upper limit in the integrals Cl has been replaced by infinity. Up 
to a normalization factor, the contribution of each k value is a moment of the potential. As it is well known [l2j, the 
leading term is vq/ p 3 for i x = l' x = 0. Expansion d is carried out for the three nucleus-nucleus potentials with 
additional Raynal-Revai transformations for the second and third terms. Analytic expansions of potentials m are 
finally obtained with 




E Vk 
n 2k 


qI a ;+ J x +3 Z -/ q* 

r k—0 r 


(16) 


where coefficients Vk are obtained from Vk after Raynal-Revai and spin coupling transformations. 

Let us evaluate coefficients Vk for 6 He=a+n+n, with the a — n potential taken from Kanada et a 1. |0. Coefficients 
vq to V 4 are given in Table [I] for J = 0 + . We also provide the amplitude of the centrifugal term 

Vcent = ^—{K + 3/2) (A" + 5/2), (17) 

zrajv 

which depends on p as 1/p 2 . It is clear from Table [I] that coefficients Vk are large and increasing with k. Integrals 
in cs must be computed with a high accuracy. Special attention must be paid to partial waves involving two-body 
forbidden states. In this case, we use a supersymmetry transform of the potential |2C|, in order to remove forbidden 
states in the three-body problem. This transformation is carried out numerically, and the resulting potential presents 
a singularity at short distances. 

From Tabled we evaluate the p value where the nuclear part is negligible with respect to the centrifugal term. In 
other words, p m ax is defined as 


N 

rmax 


= e X 


Vcent 

Pmax 


(18) 
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Values of p max are given in Tabled by assuming e = 0.01. In general they are larger for low K values for two 
reasons: (z) the centrifugal term is of course lower, and (ii) low partial waves generally involve forbidden states which 
lead to singularities in the potential, and hence to larger values of vq. 


TABLE I: Coefficients Vo to 1 T 4 in 6 He for J = 0 + , L = S = 0, and for typical partial waves (energies are expressed in MeV and 
lengths in fm). The bracketed values represent the power of 10, and 7 = l x ,£ v . 


K, 7 

K\ 7 ' 

Vo 

Vl 

V2 

V 3 

V4 

Vcent 

Pmax 

0,0,0 

0,0,0 

3.40(3) 

-7.46(3) 

-2.02(4) 

-1.53(5) 

-1.78(6) 

78 

4370 

4,0,0 

4,0,0 

1.18(3) 

-1.20(5) 

7.31(6) 

-2.13(8) 

2.87(9) 

741 

160 

8,0,0 

8,0,0 

-2.59(3) 

-1.19(5) 

5.46(7) 

-6.66(9) 

4.98(11) 

2068 

125 

4,2,2 

4,2,2 

2.61(4) 

-1.27(6) 

5.40(7) 

-1.39(9) 

1.81(10) 

741 

3520 

8,2,2 

8,2,2 

5.49(4) 

-7.82(6) 

1.06(9) 

-1.02(11) 

6.78(12) 

2068 

2660 

0,0,0 

4,0,0 

-3.41(3) 

8.04(4) 

-1.09(6) 

4.27(6) 

1.43(7) 



0,0,0 

8,0,0 

1.19(3) 

-1.08(5) 

6.21(6) 

-1.75(8) 

2.33(9) 



0,0,0 

4,2,2 

9.62(3) 

-2.41(5) 

3.47(6) 

-1.37(7) 

-4.61(7) 



0,0,0 

8,2,2 

1.40(4) 

-9.90(5) 

4.80(7) 

-1.30(9) 

1.71(10) 




From the p max values displayed in Tabled it is clear that the channel radius a of the R matrix must be very large. 
Using basis functions valid up to these distances would require tremendous basis sizes. This is solved by using a 
propagation technique which is presented in Sec. 2.3.3. 

In the analytical expansion of the potential, the maximum value k max is determined from the requirement 


'^fcmax+l „ V kmax 

a 2k max +2 ^ a 2k max • 

This yields typical values kmax ~ 3 — 4, depending on the system and on the partial wave. 


(19) 


C. Three-body R-matrix 

1. Principle of the R matrix 

The R-matrix theory is well known for many years E2 It allows matching a variational function over a finite 
interval with the correct asymptotic solutions of the Schrodinger equation. We summarize here the main ingredients 
of the R-matrix theory and emphasize its three-body aspects. The R-matrix method is based on the assumption that 
the configuration space can be divided into two regions: an internal region, with radius a, where the solution of © 
is given by some variational expansion, and an external region where the exact solutions of m are known. This is 
formulated as 


N 

XjK,int(p) = C TKi u i(p)’ ( 20 ) 

i -1 

where the N functions Ui(p) represent the variational basis, and are the corresponding coefficients. In the external 
region, it is assumed that only the Coulomb and centrifugal potentials do not vanish; we have, for an entrance channel 
l'K', 


( 21 ) 


where the amplitude is chosen as 


Ajk = i K+1 (2n/k) 5 ^ 2 , 


( 22 ) 
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and where U Jv is the collision matrix, and k = sJ^m^E/h 2 is the wave number 0 If the three particles do not 
interact, Eq. m is a partial wave of a 6 -dimension plane wave l3 


exp [i(k x .x + ky.y)] = 


M! 

(k P ) 2 


E i K J K +2{kp) 


e^i v LM L K 




For charged systems, we have 


H jk( x ) = G K+ ^(V'yK,x) ±iF K+ s(ri 7K ,x), 


(23) 


(24) 


where G^+ 3/2 and E K+3 / 2 are the irregular and regular Coulomb functions, respectively |2l|. The Sommerfeld 
parameters p n K are given by 


_ jtt m Ne 


h 2 k ’ 


(25) 


where z Jn is the effective-charge matrix HU; r) therefore depends on the channel. Notice that we neglect non-diagonal 
terms of the Coulomb potential. This is in general a good approximation as these terms are significantly smaller than 
diagonal terms 0 For neutral systems, the ingoing and outgoing functions H± K (x) do not depend on 7 and are 
defined as 


1 / 7FT \ 1/2 

H±k i x ) = ±M it) [ Jk+2(x ) ± iY K +2 {x)] , 


12 / 


(26) 


where J n (x ) and Y n (x) are Bessel functions of first, and second kind, respectively. The phase is chosen to recover the 
plane wave in absence of interaction ( U=I ). 

For bound states (E < 0), the external wave function is written as 

XjK,extiP) = W- VyK .K+2{2K,p), (27) 

where W a b{x) is a Whittaker function, and the amplitude ( k 2 = — 2itinE/ h 2 ). For neutral systems, we have 

XjK,extiP) = c Jk {kp) 1 /2 K k+2 (kp) : (28) 

where K n (x) is a modified Bessel function. 


2. The Bloch-Schrodinger equation 


The basic idea of the i?-matrix theory is to solve Eq. m over the internal region. To restore the hermiticity of the 
kinetic energy, one solves the Bloch-Schrodinger equation 


(H + C(L) - E) <F JM7r = £(L)'F JM7r 
with the Bloch operator C{L) defined as 

h 2 ^ i~« 7 A// - \ 1 / 0 L-yK 


w = 2^7 E > s (p - «°>W 


7 K 


p5/2 yg p p 


P 5/2 < y: 


7 K |i 


(29) 


(30) 


where L is a set of arbitrary constants L 7 ^-. In the following, we assume = 0 for positive energies. Formulas 
presented in this subsection are given for any channel radius a o, which can be different from a, defined in 2.3.1. 

Let us define matrix C ,/7r as 


sxJtt 

K'i' 


=< uiy'% |H + C(L) - E\ui>y^, >/, 


(31) 


where subscript / means that the matrix element is evaluated in the internal region only, i.e. for p < a q. Using the 
partial-wave expansion 0 and the continuity of the wave function at p — ao, we obtain the .R-matrix at a o from 


^K,-/'K'i a o) — 


h 2 


2?7TjvQo 


E 


*i(ao) (C J ") 


-1 


K'i' 


M' (ao). 


(32) 
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3. R-matrix propagation and collision matrix 

As shown in Sect. 2.2, the nuclear potential extends to very large distances, even for short-range nucleus-nucleus 
interactions. In other words, the asymptotic behaviour m is not accurate below distances which may be as large as 
1000 fm or more. This is a drawback of the hyperspherical method, where even for large p values, two particles can 
still be close to each other and contribute to the three-body matrix elements. 

It is clear that using basis functions valid up to distances of 1000 fm is not realistic, as the size of the basis would 
be huge. On the other hand, using a low channel radius (typically 30 ~ 40 fm) would keep the basis size in reasonable 
limits, but would not satisfy the key point of the A-matrix theory, namely that the wave function has reached its 
asymptotic behaviour at the channel radius ao- This problem can be solved with propagation techniques, well known 
in atomic physics !13j. The idea is to use ao as a starting point for the R matrix; its value is small enough to allow 
reasonable basis sizes. The A matrix is then propagated from ao to a, where the Coulomb asymptotic behaviour m 
is valid. Between ao and a, the wave functions X Jn {p) are still given by Eq. itStl . but with the potential replaced by 
its (analytical) asymptotic expansion. 

More precisely, the internal wave functions in the different intervals are given by 

N 

X^K,int(p) = C yKi u i(P) for p < a 0 , 

i= 1 

= X 7 k(p) for a 0 <p<a, (33) 

where XiK {p) are solutions of Eq. @ with the analytical expansion m of the potential term. 

The R matrix is first computed at ao with Eq. 13211 (typical values are ao ~ 20 — 40 fm). Then we consider No sets 
of initial conditions for x(p), where N 0 is the number of 7 K values (from now on we drop the Jtt index for clarity). 
We combine these sets as matrix Xo(p), and choose 

Xo(a 0 )=J, (34) 

where I is the unit matrix. 

According to the definition of the A matrix EL we immediately find the derivative at ao 

Xo'(ao) = —R 1 (ao)xo(ao) = —R 1 («o)- (35) 

a 0 a 0 

Knowing functions Xo-yK and their derivatives at ao, they are then propagated until a by using the Numerov 
algorithm |22l| . well adapted to the Schrodinger equation. The analytical form of the potential is used, with a 
summation limited to a few k values. The A matrix at a is then determined by using Eq. 13511 with Xo( a ) and x 0 (a). 
We have 

R(a) = ~Xo(a) (xo'(a)) -1 . (36) 

a 

Notice that the propagated R matrix m does not depend on the choice of Xo(ao)- In Ref- a, the propagation is 
performed through the Green function defined in the intermediate region, and expanded over a basis. The method 
presented here uses the Numerov algorithm, and does not rely on the choice of a basis. The analytical form of the 
potential in this region makes calculations fast and accurate. 

Finally the collision matrix is obtained from the R matrix at the channel radius a with 

U J r = 1 (37) 

and 

^^k,yk' = H lK (ka)8 11 ’5KK' — ka(H^ K (ka))' RJfc^ic'ia), (38) 

where the derivation is performed with respect to ka. 

Lower values of the channel radius a can be used by employing the Gailitis method |23j. In this method the 
asymptotic forms (E3) are generalized with the aim of using them at shorter distances. This means that the propagation 
should be performed in a more limited range (typical values for a are a ~ 200 — 400 fm). However this does not avoid 
propagation which, in any case, is very fast. In addition, the Gailitis method cannot be applied to charged systems, 
as it assumes from the very beginning that the coupling potentials decrease faster than 1/p. 

The extension of the .R-matrix formalism to bound states is well known for two-body systems |24| . Basically, 
the Ljk constants are defined so as to cancel the r.h.s. of Eq. <12911 . Then, the problem is reduced to a matrix 
diagonalization with iteration on the energy msi. 


4- Wave functions 


Once the collision matrix is known, the internal wave function can be determined in both intervals. Although 
the choice of Xo(ao) is arbitrary, functions x(p) entering Eq. (1331) do not depend on that choice. In the intermediate 
region Oq < p < a, functions x{p) and x 0 (p) are related to each other by a linear transformation 


X(p) = Ao {p)M. 

Matrix M is deduced by using the asymptotic behaviour m at p = a, 

X(a) = X 0 (a)M = * ext (a), 


( 39 ) 


( 40 ) 


where Xext ( c 0 is the matrix involving all entrance channels [see Eq. 4Z91- it 
Coefficients defining the internal wave function in the interval p < ao 


depends on the collision matrix, 
are finally obtained by 


„ J n — 
C~^Ki 


h XT 

2m N ^ { >yKi,YK'i' 
r y' K'i' 


d Xy K ' \ 
dp ) 

/ p=a 0 


uv{ao). 


( 41 ) 


5. The Lagrange-mesh method 


Up to now, the basis functions Ui(p) are not specified. We use here the Lagrange-mesh method which has been 
proved to be quite efficient in two-body 0 and three-body j(3] systems. Notice however that its application to 
three-body continuum states is new. 

When dealing with a finite interval, the N basis functions Ui(p) are defined as 0 


Ui(P ) 





1/2 pP N {2p/a 0 - 1) 
P 


( 42 ) 


where the Xi are the zeros of a shifted Legendre polynomial given by 


P N (2xi - 1 ) = 0 . 


( 43 ) 


The basis functions satisfy the Lagrange condition 

Ui(a 0 Xj) = (a 0 \i)~ 1/2 5ij, ( 44 ) 

where the A, are the weights of the Gauss-Legendre quadrature corresponding to the [0,1] interval, i.e. half of the 
weights corresponding to the traditional interval [-1,1]. 

The main advantage of the Lagrange-mesh technique is to strongly simplify the calculation of matrix elements 1311) 
if the Gauss approximation is used. Matrix elements of the kinetic energy (T + C) are obtained analytically jl4j . 
Integration over p provides matrix elements of the potential by a single evaluation of the potential at p = aoXi. The 
potential matrix is diagonal with respect to i and i'. 

In Ref. ;d], we applied the Lagrange-mesh technique to bound states of three-body systems. As the natural interval 
ranges from zero to infinity, we used a Laguerre mesh. It was shown that the Gauss quadrature is quite accurate for 
the matrix elements, and that computer times can be strongly reduced. 


III. APPLICATIONS 
A. Conditions of the calculations 

Here we apply the method to the 6 He and 14 Be nuclei. The a-n and 12 Be-n interactions are chosen as local 
potentials. They contain Pauli forbidden states (one in £ = 0 for a-n, and one in l = 0,1 for 12 Be-n) which should 
be removed for a correct description of three-body states For bound states, two methods are available: the 

use of a projector |27|, and a supersymmetric transformation of the nucleus-nucleus potential ]2Q]. Although both 
approaches provide different wave functions, spectroscopic properties are similar ;6|. For unbound systems, it turns 
out that the projector technique is quite difficult to apply with a good accuracy. Expansions similar to Eq. Ill (il l for 
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the projection operator provide non-local potentials. Consequently, all applications presented here are obtained with 
supersymmetric partners of the nucleus-nucleus potentials. 

As collision matrices can be quite large, it is impossible to analyze all elements. To show the essential information 
derived from the collision matrix, we rather present some eigenphases. Those presenting the largest variation in the 
considered energy range are shown. 

Analyzing the collision matrix in terms of eigenphases raises two problems. First, it is in general not obvious to link 
the eigenphases at different energies. As eigenphases cannot be associated with given quantum numbers, there is no 
direct way to draw continuous eigenphases. The procedure can be strongly improved by analyzing the eigenvectors. 
Starting from a given energy, eigenphases for the next energy are chosen by minimizing the differences between the 
corresponding eigenvectors. 

A second problem associated with eigenphases arises from the Coulomb interaction. As matrix elements of the 
Coulomb force are not diagonal, the corresponding phase shifts do not appear in a simple way, as in two-body 
collisions. Consequently, in order to extract the nuclear contribution Un from the total collision matrix U, we perform 
two calculations: a full calculation providing U, and a calculation without the nuclear contribution, providing the 
Coulomb collision matrix Uc- Then we define the nuclear collision matrix Un by 

U = U c 1 /2 U n Uc 1/2 . (45) 

As U and Uc are symmetric and unitary, the same properties hold for Un. Examples of Coulomb phase shifts will 
be given in the next sections. 


B. Application to 6 He and 6 Be 

The conditions of the calculation are those of Ref. Q. The a-n potential V a — n has been derived by Kanada 
et a 1. [r|. It contains spin-orbit and parity terms. The n-n potential is the Minnesota interaction j2f|. As bare 
nucleus-nucleus potentials cannot be expected to reproduce the 6 He ground-state energy with a high accuracy, we 
renormalize V a - n by a factor A = 1.051 (note that this value was misprinted in Ref. j6j)- This value reproduces the 
6 He experimental energy —0.97 MeV and provides 2.44 fm for the r.m.s. radius. The convergence with respect to 
K max and to the Lagrange-mesh parameters has already been discussed in Ref. 0. 

Let us first illustrate the importance of the propagation technique. In Fig. |3 we plot some elements of the J = 0 + 
collision matrix under different conditions. In each case, we compare the phase shifts for two channel radii: ao = 20 
fm and ao = 30 fm. The calculation is performed with and without propagation. For K = 0, reasonable values can be 
obtained without propagation. However, for larger K values ( I\ = 8 is displayed with i x = £ y = 0 and i x = i y = 4), 
the channel radius should be quite large to reach convergence. To keep the same accuracy, the number of basis 
functions should be increased. However, one basis function per fm is a good estimate, and this leads to unrealistically 
large basis sizes. This convergence problem is due to the long range of the potential. The propagation technique 
(performed here up to a = 250 fm) allows us to get a very high stability (better than 0.1° at all energies) even for 
rather small channel radii. Consequently calculations with high K values are still feasible. 

To illustrate the diagonalization of the collision matrix, we compare in Fig. [3 the diagonal phase shifts with the 
corresponding eigenphases. We have selected a particular case, with J = 2 + , and K max = 2. With these conditions 
the collision matrix is 4 x 4, and presents a narrow resonance near 2 MeV. In the upper part of Fig. [3 we plot 
the diagonal phase shifts. One of them presents a 180° jump, characteristical of narrow resonances. This resonant 
behaviour is also observable in two other partial waves. After diagonalization of the collision matrix (Fig. [3 lower 
part) the resonant behaviour shows up in one eigenphase only. The three other eigenphases smoothly depend on 
energy. 

The convergence with respect to K max is illustrated in Fig. [3 with the J = 0 + eigenphases. It turns out that, at 
low energies, high hypermomenta are necessary to achieve a precise convergence. However, above 4 MeV, K max = 20 
is sufficient to obtain an accuracy of 2°. 

Figure 01 gives the eigenphases for J = 0 + , 1 _ , 2 + in 6 He and 6 Be ( K max is taken as 24, 19 and 16, respectively). As 
expected, the 2 + phase shift of 6 He presents a narrow resonance. The theoretical energy (about 0.2 MeV) is however 
underestimated as the experimental value |29j | is E = 0.82 MeV. In order to provide meaningful properties for this 
state, we have readjusted the scaling factor to A = 1.020, which provides the correct energy. The 0 + and 1“ phase 
shifts show broad structures near 1.5 MeV. Similar phase shifts have been obtained by Danilin et a 1. f30l|3l| and by 
Thompson et a 1. ia with other potentials. 

In 6 Be, the ground state is found at E = 1.26 MeV with a width F = 65 keV. These values are in reasonable 
agreement with experiment [29| {E = 1.37 MeV, T = 92 ± 6 keV), the width being underestimated by the model due 
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E (MeV) 


FIG. 1: a+n+n phase shifts (J = 0 + ) for channel radii a o = 20 fm (N = 20) and ao = 30 fm (N = 30), and for different partial 
waves. Solid lines are obtained with propagation up to a = 250 fm of the R matrix (curves corresponding to different ao are 
undistinguishable), and dashed lines without propagation. 


to the lower energy. Experimentally, a 2 + state is known near E = 3.0 MeV with a width of F = 1.16 ± 0.06 MeV. 
These properties are consistent with the theoretical 2 + eigenphase, which presents a broad structure near E « 4 MeV. 
The largest Coulomb eigenphases (J = 0 + ) are shown as dotted lines in Fig.0] As expected, the Coulomb interaction 
plays a dominant role at low energies, but it cannot be completely neglected even near 10 MeV. Coulomb eigenphases 
for other spin values are very similar and therefore are not presented. Energies and widths are given in Table un 
In Table im we also present the E2 transition probability in 6 He. For the narrow 2 + resonance, we use the bound- 
state approximation. Without effective charge, the B(E 2) value for the 0 + —> 2 + transition is underestimated with 
respect to the experimental value @]. However, the E2 matrix element is very sensitive to the effective charge. A 
small correction (Se = 0.05e) provides a B(E 2) within the experimental error bars. 


C. Application to 14 Be 

As shown in previous works [fill libl Isfij. a 12 Be+n+n three-body model can provide a realistic description of 14 Be. 
The spectroscopy of the 14 Be ground state has already been investigated in non-microscopic |3J, [H, [3(| j32| and 
microscopic [38| models. Here we extend three-body descriptions to 14 Be excited and continuum states. 

The 13 Be ground state is expected to be a virtual s wave, with a large and negative scattering length (a s < —10 fm) 
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FIG. 2: Diagonal phase shifts (upper panel) and eigenphases (lower panel) for the a+n+n system (J = 2 + ,K rnax = 2). 



FIG. 3: Energy dependence of a+n+n eigenphases (J = 0 + ) for different Km ax values. 


HU- In addition, the existence of a 5/2 + d-state near 2 MeV is well established. These properties can be reproduced 
by a 12 Be-n potential 


V(r ) 


Vo + Vs i • a 
1 + exp((r - ro)/a) ’ 


(46) 


where £ is the relative angular momentum and s the neutron spin. In Eq. r 0 = 2.908 fm, a = 0.67 fm, V 0 = 43 
MeV, V s = 6 MeV. The range and diffuseness of the Woods-Saxon potential are taken from Ref. 36). The amplitudes 
Vq and V s provide E( 5/2 + ) = 2.1 MeV, and a s = —47 fm, which are consistent with the data. For the n-n potential, 
we use the Minnesota interaction, as for the 6 He study. 


With these potentials, the 14 Be ground state is found at E 


0.16 MeV, which represents an underbinding with 
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FIG. 4: Eigenphases of 6 He and e Be for different J values (solid lines). For ®Be, dotted lines represent the largest Coulomb 
eigenphases for J = 0 + . 


TABLE II: 6 He and 6 Be properties. Unless specified, experimental data are taken from Ref. [29]. 



6 He 



6 Be 


present 

exp. 

present 

exp. 

E(0+) (MeV) 

-0.97 

- 0.97 

1.26 

1.37 

T(0 + ) (keV) 



65 

92 ±6 

E( 2+) (MeV) 

0.8 

0.82 

« 4.0 

3.04 


F(2+) (MeV) 

0.04 

0.113 ± 0.020 « 1.0 1.16 ±0.06 a) Ref. tfl, b> Ref. [22] c> Ref. d) Ref. 0 

V< r 2 > (fm) 

2.44 

2.33 ± 0.04 a) 



2.57 ± 0.10 6) 



2.45 ± 0.10 c) 

B(E2,0+ -> 2+) (e 2 .fm 4 ) 

1.23 (5e = 0) 

3.2 ± 0.6 d) 


2.69 (6e = 0.05e) 


respect to experiment (—1.34± 0.11 MeV [Uj). This calculation has been performed with K max = 24, which ensures 
the convergence. The underbinding problem is common to all three-body approaches, and can be solved in two ways. 

(*) A renormalization factor A = 1.08 provides a ground-state energy at —1.30 MeV, i.e. within the experimental 
uncertainties. This procedure leads to a slightly bound 13 Be ground state, which might influence the 14 Be properties. 
(ii) A three-body phenomenological term V* *- 123 ), taken as in Ref. .ljj], i.e., 

© f (/>) = Skk' 5 W V 3 /[l + (p/p 3 ) 2 ], (47) 

reproduces the experimental energy with an amplitude V3 = 4.7 MeV (according to Ref. [fo|, we take pz = 5 fm). 
This potential is diagonal in (K, 7), and is simply added to the two-body term [see Eqs. ltTHll . ltTT)l ]. In 6 He, it was 
shown that both readjustments of the interaction provide similar results f0j|. However the renormalization factor is 
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FIG. 5: Energy of the 14 Be 0 + (squares) and 2 + (circles) states as a function of Km ax . Full symbols correspond to a renormalized 
12 Be-n potential, and open symbols to a phenomenological three-body term (see text). 


larger for 14 Be, and both methods will be considered in the following. 

The convergence with respect to K max is illustrated in Fig.E For J = 0 + , the calculations have been done with 
Kmax up to 24. The energies obtained with renormalization or with the three-body potential are very similar. This 
confirms the conclusion drawn for the 6 He nucleus §3]. 


TABLE III: Properties of the 14 Be 0 + and 2 + states. A is the renormalization factor of the 12 Be-n potential and V 3 is this 
amplitude of the three-body potential. 




A = 1.08, V 3 = 0 

A = 1 , 1/3 = 4.7 

0 + 

E (MeV) 

-1.34 

-1.34 


V< r 2 > (fm) 

3.10 

3.14 


Ps=l 

0.046 

0.033 

2 + 

E (MeV) 

-0.15 

-0.03 


V < r 2 > (fm) 

2.99 

3.04 


Ps=l 

0.192 

0.165 

0 + -*■ 2 + 

B(E2) (e 2 .fm 4 ) 

0.48 ( 5e = 0) 

0.64 (5e = 0) 



3.18 ( Se = 0.05e) 4.05 (Se = 0.05e) 


Spectroscopic properties of 14 Be are given in Tables EH and Ep The r.m.s. radii have been determined with 
2.57 fm as 12 Be radius. For the ground state, we have V< r 2 > = 3.10 frn or 3.14 fm, in nice agreement with 
experiment (3.16±0.38 fm, see Ref. 0). In all cases, the S = 1 component (denoted as Ps=i) is small (< 5%). The 
decomposition in shell-model orbitals (see Table ITYll shows that the 0 + state is essentially (« 70%) ( 2 si/ 2 ) 2 , with 
small ( 2 d 3 / 2 ) 2 and ( 2 d 5 / 2 ) 2 admixtures. 


Regarding J = 2 + , we have considered values up to K max = 16, where the number of partial waves is 172. Going 
beyond K max = 16 would require too large computer memories. Fig. E shows the energy convergence with respect 
to K max . For both potentials, the energy is below threshold, and the r.m.s. radius is close to 3 fm. A partial-wave 
analysis provides 19% of S = 1 admixture, a value much larger than in the ground state. Table HV1 suggests that the 
structure of the 2 + state is spread over many components. The ( Si/ 2 d 5 / 2 ) component is dominant (« 23%) but other 
( sd ) and ( pf) orbitals also play a role. 

E2 transition probabilities are also given in Table IIIII Without effective charge, we have B(E2,0 + —> 2 + ) = 0.48 
and 0.64 e 2 fm 4 , which is lower than for the corresponding transition in 6 He. However, the amplitudes of the proton 
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TABLE IV: Components (in %) in 14 Be wave functions. 


0+(A = 

1.08, V 3 

= 0) 0 + (A = 1, V 3 

= 4.7) 

(P 3 / 2) 2 

2.0 

2.2 


(Pi/ 2) 2 

1.0 

1.1 


Oh/ 2) 2 

70.4 

73.1 


( rf 5/2) 2 

14.6 

13.0 


(<^3/2) 2 

11.2 

9.8 


2+(A = 

1.08, V 3 

= 0) 2 + (A = 1, V 3 

= 4.7) 

P3/2P3/2 

7.7 

8.2 


P 3 / 2 / 7/2 

9.6 

9.4 


P1/2P3/2 

18.0 

18.6 


Pl/ 2 / 5/2 

5.8 

5.7 


^ 1 / 2^572 

23.2 

23.5 


^ 1 / 2^372 

19.2 

18.5 


^ 5 / 2^572 

5.6 

5.3 


^ 3 / 2^572 

4.0 

3.6 


^3/2*13/2 

3.0 

2.9 



and neutron E2 operators being even more different in 14 Be than in 6 He, the B(E2) values strongly depend on the 
effective charge. For 5e = 0.05e, we find B(E2) = 3.18 or 4.05 e 2 fm 4 according to the potential. Such transition 
probabilities should be measurable through Coulomb excitation experiments. 

In Figs. EE we present the 0 + and 2 + radial wave functions and probabilities P{x , y) defined as 

P J *(x,y) = J dCl x dQ. y x 2 y 2 | T JM7T (x, y) | 2 . (48) 



FIG. 6: Left panel: radial functions x(p) f° r the 0 + state. The curves are labeled by £ x ,£ y ,n. Right panel: probability 
P(r n n,rBe-nn), deduced from Eq. 1481 with r nn = \/2x and rs e -„„ = y/6/7y. Contour levels are plotted by steps of 0.005. 


The dominant S = 0 components are plotted. The 0 + probability shows two well distinct maxima, which resemble 
the maxima found in 6 He, corresponding to ’’dineutron” and ’’cigar” configurations. Partial waves Xjk(p) h ave 
maxima for p > 5 fm. This corresponds to distances larger than in 6 He @ where the maxima of the main components 
are located near 4 fm. As expected, the 2 + probability is similar to the 0 + probability, with two maxima. 



FIG. 7: See Fig. EJfor the 2 + state. 
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Three-body eigenphases are displayed in Fig. [8] As for the 14 Be spectroscopy the use of a three-body potential does 
not qualitatively change the phase shifts. The 1“ phase shift presents two jumps but they cannot be directly assigned 
to physical resonances. On the contrary, the 2 + phase shift shows a narrow resonance near 2 MeV. For the sake of 
completeness, 12 0+p+p mirror phase shifts are also shown in Fig. [S] As expected, no narrow structure is found. A 
very broad 0 + resonance shows up near 8 MeV, and should correspond to the 14 Ne ground state. 



E (MeV) 



E (MeV) 


FIG. 8: Three-body 12 Be-n-n and 12 0-p-p eigenphases. Solid lines correspond to a renormalized core-nucleon potential, and 
dotted lines to a phenomenological three-body term. 


IV. CONCLUSION 

In this work, we have extended the three-body formalism of Ref. Q to unbound states. As for two-body systems, 
the Lagrange-mesh technique, associated with the I?-matrix method, provides an efficient and accurate way to derive 
collision matrices and wave functions. Compared with two-body systems, three-body f?-matrix approaches are more 
tedious, owing to the coupling potentials which extend to very large distances. This behaviour is inherent to the use 
of hyperspherical coordinates which provide three-body potentials behaving as 1/p 3 , even for short-range two-body 
interactions. This problem can be efficiently solved by using propagation techniques. Here, we propagate the wave 
function and the R matrix by using the Numerov algorithm. This formalism has been extended to charged systems. 

The 6 He system has essentially been used as a test of the method, as most of its properties are available in the 
literature. The B(E 2,0 + —> 2 + ) experimental value can be reproduced with a small effective charge 6 e = 0.05e. 
We have determined a+p+p phase shifts, and found a good agreement with experiment for the 6 Be ground-state 
properties. 

Application to three-body 12 Be+n+n states is new, and has been developed in two directions. The bound-state 
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description of 14 Be provides evidence for a 2 + bound state, as expected from the shell model. The study of the 
12 Be+n+n system has been complemented by three-body phase shifts, which suggest the existence of a second narrow 
2 + resonance at E x ss 3.4 MeV. 

A limitation of the method is the slow convergence of the phase shifts with respect to the maximum hypermomentum 
Kmax■ To achieve a full convergence, values up to K rnax = 20 or more are necessary. This problem is even stronger 
for high spins, where the number of partial waves increases rapidly. A possible solution to this problem would be 
to apply the Feshbach reduction method |4Sj to scattering states. Another possible development would be to use a 
projection technique to remove Pauli forbidden states [lOj. In that case, asymptotic potentials are non local, 
which makes the calculation still heavier. 

The present model offers an efficient way to investigate bound and unbound states. In exotic nuclei, most low-lying 
states are unbound, and a rigorous analysis requires scattering conditions. The inclusion of the Coulomb interaction 
still extends the application field, and is interesting even for non-exotic nuclei. In this context, an accurate analysis 
of unbound a + a + a states seems desirable in view of its strong interest in the triple-a reaction rate ,43]. 
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